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Abstract We investigate the impact of Allee effect and dispersal on the long-term evolution of a population 
in a patchy environment, focusing on whether a population already established in one patch either successfully 
invades an adjacent empty patch or undergoes a global in-all-patch extinction. Our study is based on the 
combination of analytical and numerical results for both a deterministic two-patch model and its stochastic 
analog. The deterministic model has either two or four attractors. In the presence of weak dispersal, the 
analysis of the deterministic model shows that a high-density and a low-density populations can coexist at 
equilibrium in nearby patches, whereas the analysis of the stochastic model indicates that this equilibrium 
is metastable, thus leading after a large random time to cither an in-all-patch expansion or an in-all-patch 
extinction. Up to some critical dispersal, increasing the intensity of the interactions leads to an increase 
of both the basin of attraction of the in-all-patch extinction and the basin of attraction of the in-all-patch 
expansion. Above this threshold, while increasing the intensity of the dispersal, both deterministic and 
stochastic models predict a synchronization of the patches resulting in either a global expansion or a global 
extinction: for the deterministic model, two of the four attractors present when the dispersal is weak are lost, 
while the stochastic model no longer exhibits a metastable behavior. In the presence of strong dispersal, the 
limiting behavior is entirely determined by the value of the Allee threshold as the global population size in the 
deterministic and the stochastic two-patch models evolves as dictated by the single-patch counterparts. For 
all values of the dispersal parameter, Allee effects promote in-all-patch extinction in terms of an expansion 
of the basin of attraction of the extinction equilibrium for the deterministic model and an increase of the 
probability of extinction for the stochastic model. 

Keywords Allee effect • Deterministic model • Stochastic model ■ Extinction • Expansion • Invasion ■ 
Bistability ■ Basin of attraction • Mctastability 



1 Introduction 

Biological invasions of alien species are commonly divided into three stages: arrival, establishment, and 
expansion (Liebhold and Tobin 2008). The precise circumstances of an alien species' arrival, which refers to 
the transport of an alien species to new areas outside of its native range, are generally not known and are 
not the purpose of this article. The establishment stage refers to a growth phase of the population density 
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up to some threshold above which it is usually assumed that natural extinction is highly unlikely. However, 
if during the expansion stage, which refers to the spreading of the alien species to nearby new areas, the 
population expands in space through dispersal without significantly increasing its size, thus leading to a drop 
of its density, there might be a risk of extinction for species subject to an Allee effect. The Allee effect refers 
to a certain process that leads to decreasing net population growth with decreasing density, thus inducing 
the existence of a so-called Allee threshold below which populations are driven toward extinction. The 
causes of Allee effect identified by ecologists are numerous. They include failure to locate mates (Hopper and 
Roush 1993; Berec et al 2001), inbreeding depression (Lande 1998), failure to satiate predators (Gascoigne 
and Lipcius 2004), lack of cooperative feeding (Clark and Facth 1997), etc. Stochasticity, e.g., demographic 
and/or environmental stochasticity, may also play an important role during the critical time period when 
an alien species already established in one area starts to spread its population into a new area through 
dispersal. In this article, we think of the establishment stage as a local expansion of the population in a 
given geographical location, which involves an increase of population density in this location, while we think 
of the expansion stage as a global expansion of the population in space into nearby geographical locations 
regardless of its density. We call a global expansion successful if it leads to the population being established 
in nearby geographical locations, and unsuccessful if on the contrary the population fails to get established 
in new locations which may also lead to a global extinction (the population goes extinct in all patches) . The 
main purpose of this article is to study the critical time period when a species already established in a specific 
geographical location starts to expand in space, and determine whether the expansion stage is successful or 
not. Both Alice effect and stochasticity are central to better understand why some alien species successfully 
expand into new geographical areas, and there has been recently a growing recognition of the importance 
of these two components in biological invasions (Drake 2004; Leung et al 2004; Taylor and Hastings 2005; 
Ackleh et al 2007). Understanding their role and strength is of critical importance to gain some insight into 
why some species are more invasive than others, and may suggest some proper biological control strategies 
to regulate some populations (Liebhold and Tobin 2008). 

If an alien species subject to an Allee effect establishes its population in one area, i.e., its population is 
above the Allee threshold in this area, then the first step of population expansion is to spread to a nearby 
new area where the population is either absent or at least below the Allee threshold. A natural way to model 
this situation is to consider a two-patch model with heterogeneous initial conditions such that 

1. both patches are coupled by interacting through dispersal, and 

2. in the absence of interactions, i.e., when the patches are uncoupled, the initial conditions lead to estab- 
lishment in one patch and extinction in the other patch. 

This approach has been used previously by Alder (1993) and Kang et al (2009). In this article, we follow 
this modeling strategy to study the global expansion (population above the threshold in both patches) and 
global extinction (population below the threshold in both patches) of an alien species subject to an Allee 
effect during the critical time period between the establishment stage and the expansion stage by employing 
both a deterministic two-patch model and its stochastic analog. The objectives of our study are twofold: 
the first is to study the consequences of the inclusion of dispersal and Allee effect on the extinction and 
expansion for both deterministic and stochastic models with heterogeneous initial conditions; the second is 
to understand the effects of stochasticity by comparing the results based on both models. 

There is a copious amount of literature on the invasion and extinction of populations subject to Alice 
effects (e.g., Dennis 1989, 2002; Vcit and Lewis 1996; McCarthy 1997; Shigesada and Kawasaki 1997; Greene 
and Stamps 2001; Keitt et al 2001; Fagan et al 2002; Wang et al 2002; Liebhold and Bascompte 2003; 
Schreiber 2003; Zhou et al 2004; Petrovskii et al 2005; Taylor and Hastings 2005) which also includes various 
models in patchy environment (e.g., Amarasekare 1998a, 1998b; Gyllenberg et al 1999; Ackleh et al 2007; 
Kang et al 2009). 

In the deterministic side, Amarasekare (1998a, 1998b) investigated how an interaction between local 
density dependence, dispersal, and spatial heterogeneity influence population persistence in patchy environ- 
ments. In particular, she studied how Allee (or Allee-like) effects arise from these patchy models. Gyllenberg 
et al (1999) studied a deterministic model of a symmetric two-patch metapopulation to determine condi- 
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tions that allow the Alice effect to conserve and create spatial heterogeneities in population densities. Rather 
than exploring the global dynamics of their models, both Amarasekarc (1998a, 1998b) and Gyllenberg et 
al (1999) studied the influence of an Allee effect on local dynamics, e.g., number of equilibriums and local 
stability. There are few studies regarding the influence of an Allee effect on the extinction versus expansion 
of populations in patchy environments (e.g., Ackleh et al 2007; Kang et al 2009). Kang et al (2009) studied 
the influence of an Allee-like effect for a discrete-time two-patch model on plant-herbivore interactions where 
patches are coupled through a dispersal. Their study suggests that for a certain range of dispersal parameters 
the population of herbivores in both patches drops under the Allee threshold, thus leading to an extinction 
of the herbivores in both patches, for the majority of positive initial conditions. 

In the stochastic side, the recent work by Ackleh et al (2007) focuses on a multi-patch population model 
combining stochasticity and Alice effect. Their numerical simulations show that populations with initial sizes 
below but near their Allee threshold in each patch can still become established and invasive if stochastic 
processes affect life history parameters. The closer the population to its Alloc threshold, the greater the 
probability of invasion. A more theoretical approach based on interacting particle systems can be found in 
Krone (1999). In his model, each site of the infinite integer lattice has to be thought of as a patch which is 
cither empty, occupied by a small colony with a high risk of going extinct, or occupied by a full colony with 
a longer life span. If successful, a small colony gets established to become a full colony, while empty patches 
get colonized by a small colony due to invasions from adjacent full colonies, making space explicit. 

In this paper, although we model the population dynamics deterministically following the approach of 
Amarasekarc (1998a, 1998b), Gyllenberg et al (1999) and Ackleh et al (2007), our stochastic process as 
well as analytical results for both models are new. For the deterministic model, our focus is on the global 
dynamics of the system combining dispersal and Allee effects. In particular, we give analytical results on how 
Alice threshold and dispersal affect the geometry of the basins of attraction of the stable equilibriums. The 
stochastic model is derived from the deterministic one using a process that has two absorbing states corre- 
sponding to global extinction and global expansion, which allows to have a rigorous definition of successful 
invasion. In particular, our model is designed to study analytically the probability that a fully occupied patch 
successfully invade a nearby empty patch. To gain insight into the effects of stochasticity on the population 
dynamics, we will compare in details the results obtained for both models. 

The rest of the article is organized as follows. In section[2l we introduce the deterministic two-patch model 
with Allee effect coupled by dispersal. Based on the analysis of the invariant sets, we give a complete picture 
of the global dynamics of the system including the existence of the nontrivial locally stable equilibriums 
and the geometry of their basin of attraction. Numerical simulations have been performed to gain some 
insight into how dispersal and Allee threshold affect the exact basin of attraction of the equilibriums. In 
section O we introduce and analyze mathematically the stochastic model focusing on the time to fixation 
of the process, the existence of mctastable states and the probability of a successful invasion when starting 
from heterogeneous initial conditions. Simulations of the stochastic model have also been performed to better 
understand these aspects. In section 01 we compare the predictions based on both models, and describe the 
biological implications of our analytical and numerical results. Finally, Section [5] is devoted proofs. 



2 A simple deterministic two-patch model with Allee effects 

The first step in constructing the deterministic two-patch model is to consider single-species dynamics includ- 
ing an Allee affect as potential candidates to describe the evolution in a single patch. The two-patch model 
is then naturally derived by looking at a two-dimensional system in which both components are coupled 
through dispersal. The ecological dynamics of a single species' population subject to an Alice effect that can 
mimic the dynamics in the absence of dispersal is usually described by the model 

x = G{x)x ~ H{x) (1) 

where x(t) denotes the population density at time t. The function G measures the logistic component of 
population growth, which is given by 

G{x) = r - ax (2) 
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where r is the per capita intrinsic growth rate and a measures the extra mortality caused by intraspecific 
competition. In general, the bistability of the differential equation ([T]) is triggered by combining the negative 
density-dependence of the logistic growth G with the positive density-dependence of an additional demo- 
graphic factor represented here by the function H. The decreasing reproduction due to a shortage of mating 
encountered in low population density and the decreasing mortality due to the weakening predation risk in 
higher population density are two important examples of such factors (Stephens and Sutherland 1999) which, 
following Dercole et al (2003), can be modeled by a Holling type II functional response: H(x) = cx/(x + d). 
The resulting population model creates, under suitable parameter values, a threshold below which the pop- 
ulation goes extinct eventually and above which the population density approaches a positive equilibrium. 
The simplest and generic model that captures the population dynamics of a single species with Allee effects 
can be described by 

x = rx (x — 9)(1 — x) (3) 

where r is the per capita intrinsic growth rate after rescaling and 9 is a threshold that lies between 
and 1 after rescaling. The latter, called Allee threshold, determines whether the population goes extinct or 
establishes itself. More precisely, the population dynamics of (|3|) can be summarized as follows. 

Lemma 1 (Single species dynamics with Allee effects) If the population of a single species is described 
by ([3]), then it goes extinct when x(0) < 9 while its density goes to 1 when x(0) > 9. 

Thinking of model © as describing the population dynamics in one patch, the dynamics of two interacting 
identical patches with dispersal /i can be modeled by 

x = rx (x — 9)(1 — x) + \i (y — x) (4) 
y = ry(y-6)(l-y) + n{x-y) (5) 

where \i € [0, 1] is a dispersal parameter, representing the fraction of population migrating from one patch to 
another per unit of time. Although the system (H])-© is symmetric in x and y, asymmetry will be introduced 
by considering different initial conditions in each patch: x(0) ^ y(0). We will pay a particular attention to 
situations where one patch is initially below and the other patch above the Allee threshold, in which case, in 
the absence of dispersal, the population goes extinct in the first patch but establishes itself in the second one. 
The main objective is to understand, based on analytical and numerical results, how the dispersal parameter 
\i and the Allee threshold 9 affect the global dynamics, i.e., the limit sets of the system (HI)-© and the 
geometry of their basin of attraction. 

Our analytical results suggest the following picture of the global dynamics. Recall first that, in the absence 
of dispersal, the system has four locally stable equilibrium points, which correspond to cases when the 
population in each patch either goes extinct or gets established. In the presence of dispersal, the existence of 
(stable) limit cycles is also excluded: starting from almost every initial condition in R^_, the system converges 
to an equilibrium point. This is partly proved analytically in Theorem [1] and supported by the numerical 
simulations of Figure [TJ Therefore, we focus our attention on the existence, stability and basins of attraction 
of the equilibrium points. The effects of the dispersal parameter and the value of the Allee threshold are as 
follows. First, the dynamics of the deterministic two-patch model in the presence of weak dispersal are similar 
to that of the uncoupled system, having four locally stable equilibriums (Theorem 2]) , i.e., the extinction 
state (0,0), the expansion state (1,1) and two asymmetric interior equilibriums (x s ,y s ), (y s ,a; s ), which is 
not retained after the inclusion of demographic stochasticity, which induces two absorbing states and two 
metastable states. While increasing the dispersal parameter from 0, the basin of attraction of the extinction 
state (0,0) and expansion state (1,1) increase until a certain critical value at which both patches interact 
enough to synchronize, which drives the system to either global extinction (0,0) or global expansion (1,1): 
there are only two attractors (Theorems [5] and 2]) . Above this critical value, dispersal promotes extinction 
when the Allee threshold is below one half but promotes survival when the Allee threshold is above one half 
(Theorems [5] and [3]). Finally, in the presence of a strong dispersal, both patches synchronize fast enough so 
that the global dynamics reduce to that of a single-patch model: if the initial global density, i.e., the average 
of the densities in both patches, is below the Alice threshold then the population goes extinct whereas if it 
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exceeds the Allee threshold then the population expands globally (Theorem [3]). In other respects, for any 
value of the dispersal parameter, increasing the Alice threshold promotes extinction, and populations initially 
below the Alice threshold in both patches are doomed to extinction, whereas populations initially above the 
Allee threshold in both patches expand globally. These results are stated rigorously in the following two 
subsections. Simulation results and detailed summary are given in the last subsection. 



2.1 Global dynamics and basins of attraction 

In order to understand the global dynamics of the deterministic two-patch model, the first step is to identify 
its omega limit sets. Since the model is simply a two-dimensional ODE, its omega limit sets are either 
equilibrium points or limit cycles according to the Poincare-Bendixson theorem. As stated in the next 
theorem, when the dispersal parameter is sufficiently large, an application of the Dulac's criterion reveals 
simple dynamics by excluding the existence of limit cycles: for any initial condition, the system converges to 
an equilibrium point. 

Theorem 1 (Simple dynamics) For any c 6 [0, 3), r > and 9 € (0, 1), if 

then every trajectory of (UJ)-© converges to an equilibrium point. 

Theorem [T] indicates for instance that every trajectory of the system (jU)-© converges to an equilibrium 
point under the condition fi > r(9 2 — + l)/3 if one takes c — 0. In addition, Theorem [2] below implies that 
if limit cycles emerge for smaller values of the dispersal parameter then each of them is included in one of 
the two regions of the phase space in which the population lies above the Allee threshold in one patch but 
below the Allee threshold in the other patch, i.e., 

{(x,y) e Qq : < x < 9 and 9 < y < 1} and {(x,y) £ fl ■ < y < 9 and 9 < x < 1} 

where Qq = R^_. Numerical simulations (see Figure [IJ further suggest that, for any value of the dispersal 
parameter, there is no stable limit cycle, which implies that locally stable equilibriums are the only possible 
attractors of the system, so we focus our attention on the existence, stability and basins of attraction 
of the equilibrium points. We also would like to point out that if the system has no Alice effect, e.g., a 
metapopulation model coupled by both competition and migration with uniparental reproduction, then this 
system admits no periodic solutions (Proposition 1 in Gyllenberg et al 1999). 

It can be easily seen that the system (d))-© has three symmetric equilibriums for all positive values of the 
parameters: one boundary equilibrium given by Eq = (0,0) and two interior equilibriums given respectively 
by Eg = (9, 9) and E\ = (1, 1). For obvious reasons, we call Eq the extinction state of the system and E\ 
the expansion state. Theorem [2] below indicate that, for all parameter values, these two trivial equilibriums 
are locally stable whereas the interior equilibrium point Eg is unstable. Hence, to understand the global 
dynamics of the system, the next step is to study the geometry of the basins of attraction of the two trivial 
equilibriums, i.e., 

S = {(x(0),y(0)) e f2 ■ \im t ^oo{x{t),y{t)) = E } 
B x = {(x(0),y(0))en :Hm t ^ oo (x(t),y(t))=E 1 }. 
Letting i?o,e an d fig denote the subsets 

^0,0 = {(x, y) e f2 Q : < x < 9 and < y < 9} 
fig = {(x, y) e A) : x > 9 and y > 9} 

Lemma [1] indicates that, in the absence of dispersal, the basins of attraction of Eq and E\ for the (uncoupled) 
system are given by Bq = ft^^g \ Eg and B\ = fig \ Eg. The following theorem shows how the inclusion of a 
dispersal affects the basins of attraction. 
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Theorem 2 (Local stability and basins of attraction) 

1. The extinction state Eq and expansion state E\ are always locally stable whereas the interior fixed point 
Eg is always unstable. 

2. If 2/i > r9(l — 9) then Eg is a saddle while if 2[i < r9(\ — 6) then Eg is a source. 

3. i? o ,0 \ Eg C -Bo- If in addition 9 < 1/2 then 

B C {(x,y) en :x + y< 29}. 

4. Qg \Eg C Bi. If in addition 9 > 1/2 then 

Bi C {(x,y) G i? : x + y > 29}. 

Theorem [5] indicates that the inclusion of dispersal promotes both global extinction and global expansion 
of the system, as the basins of attraction of both equilibrium points Eq and E\ are larger in the presence 
than in the absence of dispersal. Numerical simulations further suggest that, up to a certain critical value, 
increasing the dispersal parameter translates into an increase of Bq and B\ . The value of the Alice threshold 
9 also plays an important role in the global dynamics. When the Allee threshold lies below one half, which 
is common in nature, the largest possible basin of attraction of Eq is 

{(x,y) eSl :x + y< 29}. 

Moreover, according to numerical simulations (see Figure [I}, increasing the Alice threshold promotes extinc- 
tion of the system in the sense that, the dispersal parameter being fixed, the smaller the Alice threshold, the 
smaller the basin attraction of the extinction state Eq and the larger the basin attraction of the expansion 
state Ei. Finally, we would like to point out that parts 1 and 2 of the theorem hold for the system (HJ)-© 
but not always for two-patch models with Alice effect. A counter example is provided by the metapopulation 
model coupled by both competition and migration with biparental reproduction studied by Gyllenberg et 
al (1999). The first step to prove parts 3 and 4 of Theorem [5] will be to identify the positive invariant sets 
of the system (JU)-© included into the upper right quadrant S7o = R+. Recall that a set is called positive 
invariant if any trajectory starting from this set stays in this set at all future times. Since we arc interested 
in the global dynamics of the system, our objective will be to find all the possible invariant sets in J?o- 
Notice that the union and the intersect of positive invariant sets are also positive invariant. All these positive 
invariant sets have an important role in understanding the dynamics in regions of the phase space where 
the population is below the Allee threshold in one patch but above the Allee threshold in the other patch. 
In particular, they will give us means of decomposing the phase space by restricting our attention to the 
dynamics on each invariant set and then sewing together a global solution from the invariant pieces. 

2.2 Dispersal effects and multiple attractors 

In this subsection, we study the effects of the dispersal parameter on the dynamics of the two-patch model 
when the Allee threshold is fixed. Theorem Q] suggests that the number of attractors is also equal to the 
number of locally stable equilibriums. Our study shows that the value of the dispersal parameter determines 
the number of equilibriums, thus the possible number of attractors. Let Sg denote the stable manifold of the 
unstable interior equilibrium Eg, i.e., 

Sg = {(x(0),y(0)) G Q : lim (x(t),y(t)) = Eg}. 

t— >oo 

The following theorem indicates that, when the dispersal is sufficiently large, both patches interact enough 
to synchronize, which drives the system to either global extinction or global expansion: there are only two 
stable equilibriums, the extinction state Eq and the expansion state E\. 
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Theorem 3 (Large dispersal) Assume that 



, > ^±a. P) 



Then, the system (H])-© /ias onit/ two attractors: Eq and E\. Moreover, 
1. If 9 < 1/2 and © ZioZds then 



2. If 9 > 1/2 and © ZioZds £/ien 



{(x,y) en \S e :x + y>2e} C Bj. 



{(a:,y) € f2 \S e :x + y<20} C S . 



If the inequality © holds, we can consider that the system has a very strong dispersal. Then Theorem [3] 
indicates that, when 9 < 1/2, both patches synchronize fast enough so that the global dynamics reduce to 
the one of a single-patch model: if the initial global density, i.e., the average of the densities in both patches, 
is below the Allee threshold then the population goes extinct whereas if it exceeds the Allee threshold then 
the population expands globally. In addition, the theoretical results in Theorems [2] and [3J suggest that the 
smaller the Allee threshold, the smaller the basin of attraction of the extinction state and the larger the 
basin of attraction of the expansion state. This agrees with the simulation results of Figure [T] 

Finally, in order to explore the number of locally stable equilibriums when the dispersal is small, we now 
look at the nullclines of the system. Define 

= x rx(x-9)(l-x)^ 

Then, the nullclines of the system ©-© are given by x = f(y) and y = f(x). The interior equilibriums are 
determined by the positive roots of x — f{f{x)), which is a polynomial with degree 9. This implies that the 
system has at most 8 interior equilibriums since is always a solution. 

According to the expression of the nullclines y = f(x) and x = f(y) (see Figure \S\ page [55]) , we can 
see that the number of interior equilibriums strongly depends upon the value of the dispersal parameter: 
in the presence of strong dispersal, both patches synchronize and the system has only two positive interior 
equilibriums Eg and E\ , which is confirmed by Theorem [3j while in the presence of weak dispersal, there is 
enough independence between both patches so that the system has 8 positive interior equilibriums. We are 
interested in the locally stable equilibriums since the possible number of attractors is intimately connected 
to the number of stable equilibriums. The following theorem summarizes the properties of the equilibriums 
and their stability for different parameters' values. 

Theorem 4 (Multiple attractors) 

1. If r > 0, 9,/j, G [0,1] then every trajectory converges in [0,1] 2 so all the equilibriums (x*,y*) G [0,1] 2 . 

2. If&H > r(9 2 —6 + 1) then there are only three equilibriums: Eq, Eg and E\, with Eq and E\ locally stable 
and Eg saddle. 

3. If An < r(l - 9) 2 then the nullcline 

a \ rx(x - 6>)(1 - x) 

/i 

has exactly two positive roots that we denote by < x\ < X2- Let M = maxo< a; < 2 ; 1 f(x). 

(a) If xi < M < X2 then the system has five fixed points with only two locally stable: Eq and E\. 

(b) If M > 1 then the system achieves its maximum number of equilibriums which is equal to 9; only 
four of them are locally stable: two symmetric equilibriums Eq and E\ and two asymmetric interior 
equilibriums (x s ,y s ) and (y s ,x s ). 



s 



(a) 9 = 0.30 and n = 0.10 (b) 9 = 0.40 and fi = 0.10 




Fig. 1 Solution curves and basins of attraction of the deterministic model with r = 1 and 9 < 1/2. The values of the Alice 
threshold and fraction parameter, 9 and fi, are indicated at the bottom of each simulation pictures. The dark dots are locally 
stable equilibriums. The solid lines are trajectories with arrows pointing to its converging state. The dashed line is the straight 
line: x + y = 29. The grey region is the basin attraction of the expansion state E\. The white region is the basin attraction of 
(x s , y s ) and (y s , x B ). The dark grey region is the basin attraction of the extinction state Eq. Based on Theorem fj] it is enough 
to restrict the system BJ)-© to the compact space [0, l] 2 . 
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(a) 9 = 0.70 and fi = 0.04 (b) = 0.60 and n = 0.04 




Fig. 2 Solution curves and basins of attraction of the deterministic model with r = 1 and 9 > 1/2. The values of the Alice 
threshold and fraction parameter. 8 and fi, are indicated at the bottom of each simulation pictures. The dark dots are locally 
stable equilibriums. The solid lines are trajectories with arrows pointing to its converging state. The dashed line is the straight 
line: x + y = 2d. The grey region is the basin attraction of the expansion state E%. The white region is the basin attraction of 
(x B ,y B ) and (y B ,x B ). The dark grey region is the basin attraction of the extinction state Eq. 

Part 1 of Theorem Q] suggests that we can restrict our analysis of the basin of attraction of locally stable 
equilibriums to the compact space [0,1] 2 . Moreover, from Theorem 01 we can see that when the dispersal 
parameter is small enough, the system has 9 equilibriums, including four locally stable equilibriums. In the 
presence of an Allee effect, small dispersal may promote survival: patches that are below the Alice threshold 
are rescued by immigrants from adjacent patches above the Allee threshold. This implies that when dispersal 
is introduced to a system with an Alice effect, populations can exist at intermediate densities, corresponding 
to the equilibriums (x s ,y s ) and (y s ,x s ), as a source-sink system, or expand to high density E±. Moreover, 
according to perturbation theory (Simon 1974; Amarasckarc 2000), both asymmetric interior equilibriums 
appear from the equilibriums (0,1) and (1,0) of the uncoupled system, i.e., in the absence of dispersal, 
caused by the small perturbation fi. Therefore, we have x s = 0([i) and y s = 1 — 0(/j,). Finally, note that the 
absence of limit cycles given by Theorem [T] when (J6j> holds combined with Theorem [4] implies that 

Corollary 1 (Four attractors) If the system has four locally stable equilibriums and inequality ^ 

holds for some c € [0, 3), then the system has exactly four attractors. 

When the system has four attractors as stated in Corollary [TJ the simulations in Figure [T] suggest that the 
smaller the dispersal, the smaller the basin of attraction of the extinction state and the expansion state, but 
the larger the basin attraction of the asymmetric interior equilibriums. In particular, if /j. — 0, then 

B ^n . e , B x ^n e , B s -)- R 2 + \ {Q 0i g U Q e ) 

where B s denotes the basin of attraction of the asymmetric interior equilibriums. 



2.3 Simulations and Summary 

Theorem [3] suggests that the larger the Allee threshold, the larger the basin of attraction of the extinction 
state and the smaller the basin of attraction of the expansion state when inequality (|7|) holds. The simulations 
shown in Figure [T] confirm this and give us a more complete picture of how the dispersal /i and Alice 
threshold 6 affect the exact basin of attraction of the locally stable equilibriums including asymmetric 
interior equilibriums: 
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Effects of dispersal fi - Fix Allee threshold 8 and growth rate r, let dispersal \i vary. 

1. When fj, is small so that the system has four locally stable equilibriums (ii smaller than some critical 
value [i c ) , the smaller the dispersal, the smaller the basin of attraction of the extinction state and the 
expansion state, but the larger the basin of attraction of the asymmetric interior equilibriums (see 
(d) and (f) of Figure [J). This indicates that smaller dispersals promote persistence of the populations 
in both patches by creating sink-source dynamics. 

2. When fj, is large so that the system has only two attractors Eq and E\ (fx larger than the critical 
value fx c ), the larger the dispersal, the larger the basin of attraction of the extinction state but the 
smaller the basin of attraction of the expansion state when < 1/2 (see (a) and (c) of Figure [l}. 
When 8 > 1/2, the monotonicity is flipped due to the symmetry of the system (see Figure [2]). 

3. Extreme cases: when fx is very small, the two-patch model behaves nearly like the uncoupled system, 
having four attractors and almost the same basins of attraction, while when fi/r— >oo. the global 
population behaves according to a one-patch system with Allee threshold 2(9, in particular 

B — ► {(x,y)eM 2 + :x + y<26}. 

Effects of Allee threshold 8 - Fix dispersal fi and growth rate r, let Allee threshold 8 vary. 

Regardless of the number of locally stable equilibriums, the larger the Allee threshold, the larger the 
basin of attraction of the extinction state but the smaller the basin of attraction of the expansion 
state (see (a), (b), (e) and (f) of Figure [J) . 

Effects of Growth rate r - Fix dispersal fi and Allee threshold 8, let growth rate r vary. 

By introducing the new time r = t/r, we can scale off the parameter r of the system ((U)-© so the 
dispersal fi becomes ji/r. This implies that the growth rate r and the dispersal parameter fi have 
opposite effects on the basin of attraction of the locally stable equilibriums, i.e., increasing the value 
of r is equivalent to decreasing the value of [i. 

Tables [T][2] give a complete picture, based on our analytical and numerical results, of how dispersal and Alice 
threshold affect the basin of attraction of the locally stable equilibriums. We only focus on the case 8 < 1/2 
but similar results can be deduced when 6 > 1/2 using the symmetry of the system (©. 



Table 1 Summary for the deterministic model when 9 < 1/2 and variations of the parameters are restricted to the case when 
the system has only two attractors Eg and E\ . 



Two attractors and 9 < 1/2 


Parameters 


Basin of attraction of Eq 


Basin of attraction of E\ 


Dispersal fi "J" 


Bo t 


B\ I 


Alice threshold 9 f 


Bo t 


Bi I 


In particular, if fi/r oo then Bo — > {(x,y) £ : x + y < 29}. 



Table 2 Summary for the deterministic model when 9 < 1/2 and variations of parameters are restricted to the case when the 
system has four attractors: Eq,Ei and (x a ,y B ), (y s ,x s ). 



Four attractors and 9 < 1/2 


Parameters 


Basin of attr. of Eo 


Basin of attr. of asymmetric equilibriums 


Basin of attr. of E\ 


Dispersal /i J. 


Bo I 


B s t 


Bi I 


Alice threshold 9 I 


Bo I 


no monotonicity 


Bi t 


-* 


Bo — > ^0,8 


B s -> \ (% e U fig) 


B\ — > Qg 
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3 Definition of the stochastic model and main results 

While the deterministic model is similar to the one in Ackleh et al (2007), our stochastic model differs 
from theirs, which is derived naturally from the deterministic model by including independent Poisson 
increments, i.e., variability in birth, death and migration events. This gives rise to a multi-patch individual- 
based model for which they study numerically the probability of a successful invasion, defined as the event 
that the population size in one patch exceeds some denominated threshold. However, well-known results 
about irreducible Markov chains imply that the population is driven almost surely to extinction which 
corresponds to the unique absorbing state of their stochastic process. In contrast, we model stochastically 
the two-patch system via a process that has two absorbing states corresponding to a global extinction and a 
global expansion, respectively. This allows to have a definition of successful invasion more rigorous and more 
tractable mathematically. In particular, while their stochastic model is designed to study numerically the 
probability that a population starting near the Allee threshold in each patch gets successfully established, 
our model is designed to study analytically the probability that a fully occupied patch successfully invade 
a nearby empty patch. More precisely, to understand the effect of stochasticity on the interactions between 
both patches, we introduce a Markov jump process that, similarly to the deterministic model, keeps track 
of the evolution of the population size in each patch. To obtain a Markov process, the state is updated 
at random times represented by the points of a Poisson process with a certain intensity making the times 
between consecutive updates independent exponentially distributed random variables. Motivated by the fact 
that the unit square S = [0, l] 2 is positive invariant for the deterministic model, we will choose this set 
as the state space, i.e., the state at time t is a random vector r\ t = (Xt,Y t ) £ S, where the first and 
second coordinates represent the population size in the first and second patch, respectively. Following the 
deterministic model, the stochastic dynamics involve three mechanisms: expansion, extinction, and migration. 
To model the presence of an Allee affect, we again introduce a threshold parameter G (0, 1) that can be seen 
as a critical size under which the population undergoes extinction and above which the population undergoes 
expansion, i.e., Allee threshold. This aspect is modeled by assuming that each component of the stochastic 
process jumps independently at rate r > to either (extinction) or 1 (expansion) depending on whether 
it lies below or above the Allee threshold. Recall that an event "happens at rate r" if the probability that 
it happens during a short time interval of length At approaches r{At) as At — ¥ 0. In particular, expansion 
and extinction are formally described by the conditional probabilities 

P (X t+At = l\X t >6) = P (Y t+At = l\ Y t >0) = rAt + o(At) 
P (X t+At = 0\ X t <6) = P (Y t+At = | Y t < 6) = rAt + o(At). 

This is also equivalent to saying that the waiting time for an expansion or an extinction is exponentially 
distributed with mean 1/r. Given that the population size in a given patch is at the Alice threshold, we 
flip a fair coin to decide whether an expansion or an extension event occurs at that patch which, in view of 
well-known properties of Poisson processes, implies that 

P (X t+A t = l\X t = 6) = P {Yt+At = l\ Y t = 6) = (r/2)At + o(At) 
P {Xt+At = 0\X t = 6) = P {Yt+At =Q\Y t = 0) = (r/2)At + o(At). 

To understand the effects of inter-patch interactions on the evolution of the system, we also include migration 
events consisting of the displacement of a fraction fi of the population of each patch to the other patch. We 
assume that these events occur at the normalized rate 1, therefore migrations arc described by 

P{{X t +At,Yt+At) = (l-»){Xt,Yt)+ti{Yt,Xt)) = At + o{At), 

We refer to Figure[3]for a schematic illustration of the dynamics, where dark rectangles represent parts of the 
populations which are interchanged in the event of a migration. To analyze mathematically the stochastic 
process, it will be useful to look at the model as a simple example of interacting particle system. Interacting 
particle systems are continuous-time Markov processes whose state space maps the vertex set of a connected 
graph into a set representing the possible states at each vertex. The evolution is described by local interactions 
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Fig. 3 Schematic representation of the stochastic model r)t = (Xt,Yt). The dark rectangles represent parts of the populations 
which are exchanged in the event of a migration. 



as the rate of change at a given vertex only depends on the configuration in its neighborhood. In particular, 
the Markov process {r]t}t can be seen as an interacting particle system evolving on a very simple graph 
that consists of only two vertices, representing both patches, connected by one edge, indicating that patches 
interact. The reason for looking at the stochastic model as an example of interacting particle system is that 
this will allow us to construct the process graphically from a collection of independent Poisson processes 
based on an idea of Harris (1972), which is a powerful tool to analyze the process mathematically. 

We now describe in details the behavior of the process along with our main results. Note that, considering 
a stochastic model rather than a deterministic one, the long-term behavior is described by a set of invariant 
measures on the state space rather than single point equilibriums. To the two trivial equilibriums of the 
deterministic model, Eq and E\ , correspond two invariant measures which are Dirac measures that concen- 
trate on those two points, respectively. These two measures are two absorbing states: the configuration in 
which both patches are empty and the configuration in which both patches are fully occupied. We call global 
extinction and global expansion the events that the process eventually fixates to the first and the second 
absorbing state, respectively. Interestingly, to the two asymmetric equilibriums of the deterministic model 
in the presence of weak dispersal correspond two quasi-stationary distributions representing two metastable 
states of the stochastic process (see Theorem [7]): depending on the initial configuration, the transient behav- 
ior might be described by one of these two quasi-stationary distributions, but after a long random time in the 
presence of weak dispersal (see Theorem[6|), the system fixates to one of the two absorbing states, suggesting 
that situations predicted by the deterministic model in which a small population can live next to a large 
population are artificially stable. Another important question is how stochasticity affects the geometry of the 
basins of attraction of the two absorbing states, although strictly speaking there is no basin of attraction for 
the stochastic model since the limiting behavior might be unpredictable, and how fast the system fixates. We 
will see that there is a set of initial configurations for which the limiting behavior of the stochastic process 
is predictable, and fixation to one of the two absorbing states occurs quickly (see Theorem [5]) • Starting from 
any other configuration, the limiting behavior becomes unpredictable in the sense that the process may reach 
any of the two absorbing states with positive probability. In the presence of weak dispersal, however, the 
limit is almost predictable in the sense that the probability that the system undergoes a global expansion 
after exiting one of its metastable states approaches zero or one (see Theorem [5|). Whether the system fixates 
to one or the other absorbing state strongly depends on the value of the Allee threshold. The limit is less 
and less predictable and the time to fixation shorter and shorter as the dispersal parameter increases. 
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3.1 Predictable behavior 

In order to describe rigorously the behavior of the stochastic model introduced above, our main objective is 
to estimate the times to fixation 

t+ = inf {t > : X t = Y t = 1} and r" = inf {t > : X t = Y t = 0}, 

and the corresponding probabilities of fixation, 

P{t = t + ) and P (t = t~) where r = min(r + , r - ), 

as a function of the initial configuration and the three parameters of the system. As previously explained, 
in contrast with the deterministic model which can have up to four distinct attractors, with probability one, 
either global expansion or global extinction occurs for the stochastic process, i.e., 

P(t<oo) = P(t = t + ) + P(t = t-) = 1. 

The state space can be divided into four subsets. Starting from only two of these subsets the limit is 
predictable in the sense that 

P(r = r+) G {0,1}. 

We call an upper configuration any configuration of the system in which the population size in each patch 
exceeds the Alice threshold, and a lower configuration any configuration in which the population size in each 
patch lies below the Allee threshold. These sets are denoted respectively by 

fl+ = {(x,y) G S : x > 9 and y > 9} = fig tl \ {(x,y) : x = 9 or y = 9} 
fi~ = {(x,y) G S : x < 9 and y < 9} = Q a , g \ {(x,y) : x = 9 or y = 9}. 

Note that the set of upper configurations is closed under the dynamics, i.e., once the system hits an upper 
configuration, the configuration at any later time is also an upper configuration. This implies that, starting 
from an upper configuration, global expansion occurs with probability one. Similarly, starting from a lower 
configuration, global extinction occurs with probability one. By representing the process graphically, the 
time to fixation can be computed explicitly, as stated in the following theorem. 

Theorem 5 (time to fixation) We have 

6r + 1 

E[r+| (X ,Y ) e Q + ] = E[t~| (X ,Yo)eQ-} = — 

Zr z 

The previous theorem indicates that, starting from an upper configuration, the system converges with prob- 
ability one to the absorbing state (1,1), whereas starting from a lower configuration, it converges with 
probability one to the other absorbing state (0,0). This result can be seen as the analog of Theorem [5] which 
states that the sets of upper and lower configurations are included in the basin of attraction of the equilib- 
rium points Ei and Eq, respectively. Theorem [S] also indicates that, when the rates at which expansions, 
extinctions, and migrations occur are of the same order, the expected time to fixation is quite short. 

3.2 Metastability 

The long-term behavior of the process starting from a configuration which is neither an upper configuration 
nor a lower configuration is more difficult to study as the probabilities of global expansion and global 
extinction are both strictly positive, which we shall refer to as unpredictable behavior. We will prove that, in 
any case, the system hits either an upper or a lower configuration at a random time which is almost surely 
finite, after which it evolves as indicated by Theorem [S] Hence, the time to fixation and probabilities of 
global expansion and extinction can be determined by estimating the hitting times 

T+ = inf {t > : (X u Y t ) £ Q + ) and T~ = inf {t > : (X t , Y t ) G J?~} 
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and the corresponding hitting probabilities 

P (T = T+) and P(T = T~) where T = min(T+, T~), 
since Theorem [5] implies that 

fir 4- 1 

E[rl = E[T] + ^— V- and P(t = t+) = P(T = T+). 
2r 2 

Even though our next results hold for any values of the parameters, they indicate that interesting behaviors 
emerge when the dispersal parameter \x is small. In contrast with the deterministic model which, in this case, 
has four attractors, as indicated by Theorem 2J the stochastic model first exhibits a mctastablc behavior 
by oscillating for an arbitrarily long time around one of the two nontrivial equilibriums of the deterministic 
model, and then fixates to one of its two absorbing states. The limit is almost predictable as the probability of 
global expansion approaches either or 1 depending on the value of the threshold parameter. For simplicity 
and since the system is symmetric, we shall assume that Xq = and Yo = 1 but the proofs of our results 
easily extend to the more general case when 

< n « min{|X o -0|,|y o -0|}. 

Recall that, starting from an upper configuration or a lower configuration, the time to fixation is rather 
small. In contrast, when Xq = 0, Yo = 1 and fi is small, the stochastic process converges to a quasi- 
stationary distribution in which the population size at patch X is relatively close to and the population 
size at patch Y relatively close to 1, and stays at its quasi-stationary distribution for a very long time, i.e., 
the expected value of T is large. However, due to stochasticity, the system reaches eventually an upper or a 
lower configuration, and then fixates rapidly. The next theorem gives an explicit lower bound of the expected 
value of the hitting time, which is the time the system stays at its quasi-stationary distribution. 

Theorem 6 (metastability) For any initial configuration, we have 

P(T < oo) = P(t < oo) = 1. 

Moreover, 

E[r|(* ,y o) = (o,i)] > ^(iif)" 1 

where 

imn(rn(l-fl),ln(fl)) 

Mi - ti) _ ' 

Note that, when the Allee threshold is bounded away from and 1, and the dispersal parameter is small, 
no is large, and so is the expected value of the hitting time T. Note also that, before the hitting time, no 
expansion event can occur at patch X while no extinction event can occur at patch Y. This indicates that 
the metastable state of the stochastic two-patch model is described by the stationary distribution of the 
Markov process fj t = (X t ,Y t ) with state space S = [0, l] 2 , and whose evolution is given by 

P (X t+At = 0\X t ^0) = P (Y t+At = 1 | Y t ^ 1) = rAt + o(At) 
P((X t+At ,Y t+At ) = (X t ,Y t )+ f x(Y t ,X t ))=At + o(At), 

where At is a small time interval. That is, the process {fjt}t is obtained from {r] t }t by assuming that only 
extinction events at patch X and only expansion events at patch Y can occur, which indeed describes the 
evolution of the original process {n t }t before it reaches an upper or a lower configuration. Letting v denote 
the stationary distribution of this new process, the behavior of the stochastic two-patch model before the 
hitting time T is described by the following theorem. 
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Theorem 7 (metastable state) Under the measure v we have 



and E„(Y t ) > 1 



This indicates that, when \i is small, the population size at patch X is close to (i.e., 0(h)) and the 
population size at patch Y close to 1 (i.e., 1 — 0(h)). The expected values above have to be thought of as the 
analog of the two asymmetric equilibriums of the deterministic model: (x s ,y s ) and (y s ,x s ). After evolving 
a long time according to the quasi-stationary distribution u, the process hits either an upper or a lower 
configuration, so the last question we would like to answer is whether global expansion or global extinction 
occurs after the system exits its metastable state. Starting from an upper or a lower configuration, the answer 
is given by Theorem [SJ Starting from Xq = and Yq = 1, the symmetry of the model implies that 



P(t 



P(t 



1/2 whenever 9 = 1/2. 



Our last result shows that, when 9 ^ 1/2 and > is small, the limiting behavior of the system is almost 
predictable in the sense that the probability of global expansion approaches either or 1. 

Theorem 8 (hitting probabilities) Assume that 9 < 1/2. Then 



where 



P(T = T~ \ (X ,Yq) = (0,1)) 
P(T = T+ | (X ,F ) = (0,1)) 

Ln(20) 



< 



m 



_HI-h)_ 

The previous theorem indicates that, when /i > is small, 

P (r = r+ | (X , Y„) = (0, 1)) = P(T = T+\ (X ,Y ) = (0, 1)) 
= 1 - P(T = T~ | (X ,Y ) = (0,1)) > l-(l + r)- m ° s 



1. 



In particular, in contrast with the deterministic model for which the limit depends on the initial condition 
and the geometry of the basins of attraction, starting from any initial configuration but an upper or a lower 
configuration, the limiting behavior of the stochastic model is only sensitive to the value of the parameters, 
with the Allee threshold 9 playing a central role. 



3.3 Simulation results 

While Theorem [5] gives an exact estimate of the time to fixation starting from particular initial conditions, 
the other results provide theoretical lower and upper bounds that allows us to gain a valuable insight into the 
long-term behavior of the stochastic two-patch model in the presence of weak dispersal. To better understand 
the combined effect of the Allee threshold and dispersal parameter when starting from heterogeneous initial 
conditions, we refer the reader to the numerical simulations of Figure 0] and Tables [3H] The left panel of the 
figure represents the probability of a global extinction, with the probability increasing with the darkness, 
and the right panel the expected time to fixation, with time increasing with the darkness, as a function of 
the dispersal parameter and the Allee threshold. The tables provide some numerical values of the probability 
of extinction and expected time to fixation averaged over 10,000 independent realizations of the stochastic 
process for specific values of the parameters. The predictions based on Theorems [5] and [H] that the time 
to extinction blows up and the probability of extinction approaches either zero or one in the presence of 
weak dispersal appears clearly looking at the left side of both panels and the left column of the tables for 
which h — 0.02. The left panel and Table [3] further indicate that the probability of a global extinction 
depends non-monotonically upon the dispersal parameter: when the Allee threshold is below one half, the 
probability of extinction first increases with the dispersal parameter and then decreases after the dispersal 
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probability of extinction time to fixation 




0.02 0.1 0.2 0.3 \t 0.5 0.02 n 0.2 



Fig. 4 Simulation results for the probability of a global extinction and the time to fixation of the stochastic model starting 
with one empty patch and one fully occupied patch and with growth parameter r = 0.25. Left: the gradation of grey represents 
the probability of a global extinction ranging from = white to 1 = black. Right: the gradation of grey represents the time 
to fixation ranging from = white to 100 or more = black. In both pictures, the probability and time are computed from the 
average of 10,000 independent simulation runs for 200 different values of the Allee threshold ranging from 0.25 to 0.75. These 
are further computed for 190 different values of the dispersal parameter ranging from 0.02 to 0.50, and 76 different values of 
the dispersal parameter ranging from 0.02 to 0.20, respectively. 

reaches a critical value that depends on 9, which can be easily seen in the row 9 = 0.45 of the table. When 
the Allee threshold exceeds one half, the monotonicity is flipped. Simulations also indicate that, the dispersal 
parameter being fixed, the probability of extinction increases as the Allee threshold increases. Although we 
omit the details of the proof, this can be easily shown analytically invoking a standard coupling argument 
to compare two processes, the first one with Allee threshold 9\ and the second one with #2 > 9%, the other 
parameters being the same for both processes. The black triangle labeled 1 in the upper right corner of the 
left picture reveals that global extinction occurs almost surely when 9 > 1 — fi. Indeed, starting from the 
heterogeneous condition Xq — and Yq = 1, after the first migration event, we have 

Xt = H and Y t = 1 — fi and so max(X t ,Y t ) = 1 — fj, < 9. 

In particular, both patches are below the Allee threshold from which it follows that the population goes 
extinct eventually. Almost sure global expansion in the parameter region corresponding to the lower right 
white triangle labeled 2 can be proved similarly. Finally, as suggested by Theorem [SJ the right picture 
and Tabic |4] indicate that the expected value of the time to fixation increases as the dispersal parameter 
decreases but also as the Allee threshold gets closer to one half, which can again be proved analytically based 
on standard coupling arguments even through we omit the details of the proof. 

4 Comparison and biological implications 

Recall that, in the absence of interactions between patches, both the deterministic model and the stochastic 
model predict a local expansion in patches where the initial population size is above the Allee threshold and 
a local extinction in patches where the initial population size is below the Alloc threshold. This induces the 
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Table 3 Probability of extinction (r = 0.25) 





0.020 


0.050 


0.100 


0.200 


0.300 


0.400 


0.500 


6 = 0.75 


1.000 


0.999 


0.995 


0.999 


1.000 


1.000 


1.000 


6 = 0.70 


1.000 


0.993 


0.978 


0.991 


0.972 


1.000 


1.000 


6 = 0.65 


0.999 


0.973 


0.928 


0.928 


0.944 


1.000 


1.000 


e = o.60 


0.992 


0.916 


0.834 


0.823 


0.905 


0.852 


1.000 


e = 0.55 


0.910 


0.759 


0.681 


0.719 


0.779 


0.857 


1.000 


e = o.50 


0.439 


0.506 


0.496 


0.502 


0.502 


0.488 


0.000 


e = 0.45 


0.056 


0.247 


0.316 


0.288 


0.230 


0.141 


0.000 


e = 0.40 


0.004 


0.083 


0.163 


0.176 


0.099 


0.000 


0.000 


9 = 0.35 


0.000 


0.030 


0.066 


0.078 


0.062 


0.000 


0.000 


e = o.30 


0.000 


0.007 


0.022 


0.011 


0.000 


0.000 


0.000 


6> = 0.25 


0.000 


0.001 


0.007 


0.002 


0.000 


0.000 


0.000 



Table 4 Time to fixation (r = 0.25) 





0.020 


0.050 


0.100 


0.200 


0.300 


0.400 


0.500 


e = 0.75 


152.711 


28.888 


19.749 


16.133 


15.165 


14.994 


15.058 


e = o.70 


334.855 


38.060 


21.694 


16.526 


16.291 


15.183 


14.960 


e = 0.65 


815.171 


52.247 


24.436 


17.526 


16.608 


15.027 


15.164 


e = o.60 


2052.663 


70.766 


27.399 


18.973 


16.458 


16.173 


15.065 


e = 0.55 


6058.586 


91.957 


29.515 


19.258 


16.877 


16.092 


14.820 


e = o.50 


10520.799 


102.694 


30.245 


20.034 


18.747 


18.606 


14.884 


e = 0.45 


5566.830 


91.369 


29.453 


19.086 


16.838 


15.946 


14.907 


e = 0.40 


2075.433 


70.278 


27.318 


18.892 


16.220 


14.838 


14.946 


e = 0.35 


829.009 


51.486 


23.954 


17.500 


16.421 


14.973 


14.936 


e = o.30 


339.504 


37.717 


21.467 


16.309 


14.867 


15.143 


14.993 


e = 0.25 


149.811 


28.528 


19.424 


16.066 


15.021 


14.885 


14.934 



Table 5 Comparison between deterministic and stochastic models 





Deterministic model 


Stochastic model 


Dispersal parameter 


9 < 1/2 | 9 > 1/2 


6 < 1/2 | 6» > 1/2 


No dispersal /i = 


4 attractors 


4 absorbing states 


Weak dispersal fi > 


4 attractors 
Bq and B\ f as fi "f fi c 


2 absorbing states 
+ 2 metastable states 
P (expansion) as 1 \ P (extinction) ca 1 


Critical dispersal /i = /j, c 


Both patches synchronize 
4 attractors — > 2 attractors 


2 absorbing states 


Stronger dispersal fi > /i c 


2 attractors 
Bo f as fj, t | Bi t as M T 


2 absorbing states 
unpredictability + quick fixation 


Very strong dispersal fj,/r large 


same behavior as one-patch model 


same behavior as one-patch model 
when starting from (0, 1) 



existence of four locally stable equilibriums for the deterministic model, and four absorbing states for the 
stochastic model, which correspond to cases when the population in each patch either goes extinct or gets 
established. Including interactions between patches, our results for the deterministic model indicate that, in 
the presence of weak dispersal, the dynamics retain four attractors, just as in the absence of interactions, up 
to a critical value \i c when the patches synchronize: the two asymmetric equilibrium points are lost so that 
only global expansion and global extinction can happen. In contrast, including both stochasticity and even 
weak interactions, only the two absorbing states corresponding to global expansion and global extinction 
are retained. The most interesting behaviors emerge when the dispersal is weak, in which case, to the two 
asymmetric locally stable equilibriums of the deterministic model, correspond two metastable states for the 
stochastic model. 
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Looking at the global dynamics, the predictions based on the analysis of the deterministic two-patch 
model indicate that below the critical value fi c dispersal promotes global expansion and global extinction in 
the sense that the basins of attraction of the two trivial fixed points expands while increasing the dispersal 
parameter. Above the critical value \x c dispersal promotes a global expansion when the Allee threshold exceeds 
one half but promotes global extinction in the more realistic case when the Allee threshold lies below one 
half. As mentioned above, in the presence of weak dispersal, both asymmetric equilibrium points become 
two metastable states, i.e., quasi-stationary distributions, after the inclusion of stochasticity, suggesting 
that situations in which a small population lives next to a large population are artificially stable: in such 
a context, the two-patch system evolves first as dictated by one of the two quasi-stationary distributions 
then, after a long random time, experiences either a global expansion or a global extinction. In addition, 
the long-term behavior of the stochastic model becomes almost predictable in the sense that, with very high 
probability, the system will undergo a global expansion when the Alice threshold lies below one half and a 
global extinction when the Allee threshold exceeds one half, which is of primary importance to predict the 
destiny of heterogeneous two-patch systems in the presence of weak dispersal. While increasing the dispersal 
parameter, the stochastic model no longer exhibits a metastable behavior, the time to fixation decreases, and 
the long-term behavior becomes more and more unpredictable. In the presence of a very strong dispersal, 
however, the analysis of the deterministic model and the stochastic model starting from a heterogeneous 
configuration give the same predictions. In this case, both patches synchronize enough so that the global 
dynamics reduce to that of a single-patch model: if the initial global density, i.e., the average of the densities 
in both patches, is below the Allee threshold then the population goes extinct whereas if it exceeds the Alice 
threshold then the population expands globally. 

Our analysis of idealized two-patch models is an important first step to understand more realistic multi- 
patch systems. Empirical data indicate that Allee thresholds in nature vary accross species and habitat 
types but are typically much smaller than one half. The predictions, based on the deterministic model in 
the presence of enough dispersal so that patches synchronize and on the stochastic model in the general 
case, that populations usually expand successfully when the Allee threshold is small is due to the fact that 
only two patches interact. Literally, the critical threshold 1/2 has to be thought of as one divided by the 
number of patches. Looking at a multi-patch model in which n patches interact all together, our analytical 
results suggest that a critical behavior should emerge for Allee thresholds near l/n when starting with a 
population established in only one patch, and more generally the number of patches where the population 
is initially established divided by the number of interacting patches. Therefore, even for realistic values of 
the Alice threshold, the long-term behavior is no longer straightforward in the presence of a large number of 
patches. Numerical simulations can also provide a valuable insight into the long-term behavior of multi-patch 
models including additional refinements such as density-dependent dispersals, heterogeneous environments 
with possibly different Alice thresholds in different patches, and more importantly the inclusion of a spatial 
structure through a network of interactions represented by a two-dimensional regular lattice or more general 
planar graphs rather than a complete graph where patches interact all together. 



5 Proofs 

Preliminary results 

As previously explained, the key to proving our main results is to first identify a number of sets which are 
positive invariant for the system (jH)-©. This will they give us means of decomposing the phase space by 
restricting our attention to the dynamics on each invariant set and then sewing together a global solution from 
the invariant pieces. Our first preliminary result indicates that, starting from any biologically meaningful 
initial condition, that is any condition belonging to Q$ := the trajectory of the system stays in the 
upper right quadrant and is bounded. 

Lemma 2 The system (H])-© is positive invariant and bounded in 
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Proof Assuming by contradiction that the system (H])-© is not positive invariant in upper right quadrant, 
we can find xq , Do > and a time T > such that 

x(0) = x and 2/(0) = y implies {x(T),y(T)) £ Q . 

Let F denote the boundary of J7o, i.e., 

r = {{x, y) e M 2 : (x = and y > 0) or (x > and y = 0)}. 

By continuity of the trajectories, the intermediate value theorem implies the existence of a time t < T such 
that (x(t),y(t)) G r therefore 

S := sup{t<T:(x(t),y(t))er} 
is well defined and (x(t),y(t)) ^ J?o for all t £ (S, T}. Then, we have the following alternative. 

1. If x(S) = y(S) = then x{t) = y(t) = for all t > S, which contradicts the existence of T. 

2. If x(S) = and y(S) > then x'(S) = fiy(S) > so 

there exists e > such that x(t) > for all t £ (S,S + e). 

This contradicts the existence of S. 

3. If x(S) > and y(S) = 0, the same argument exchanging the roles of the functions x and y leads again 
to a contradiction. 

In conclusion, if x(0) > and y(0) > then (x(t),y(t)) G £2o at all positive times t, which establishes the 
first part of the lemma. This also implies that, starting from any initial condition in i?0i 

x + y = rx (x — 9)(1 — x) + ry (y — 8)(1 — y) 

= r [-(x 3 + y 3 ) + (l + 6)(x 2 +y 2 )-6{x + y)} < 

whenever x + y is larger than some M(8) > 0. Therefore, 

max(x(i), y(i)) < x(t) + y(t) < M(9) for alH large enough. 

This completes the proof of Lemma [21 I 

It follows from the previous lemma that, excluding the initial condition in which both patches are initially 
empty, the population densities in both patches are simultaneously positive at any positive time. This implies 
in particular that the trivial equilibrium Eq is the only boundary equilibrium. 

Lemma 3 // (x(0), 2/(0)) eR 2 + \ {(0, 0)} then x(t) > and y(t) > for all t > 0. 

Proof By symmetry, we may assume that x(0) > and y(0) > 0. We first apply Lemma [5] to get 

< x(t),y(t) < M := max(M(#), x(0), y(0)) at all times t > 

where M(9) is as in the proof of Lemma [5J In particular, 

x'(t) = rx(x-6)(l-x) + /j,(y-x) > [r(x - 6){1 - x) - fj\ x 
> - [r max (6,(M-6)(M - 1)) + n] x > -Kx 
for some constant K < oo. Therefore, 

x(t) > x(0) exp(-Kt) > for all t > 0. 
Finally, if j/(0) > then the same holds for y(t), while if j/(0) = then 

y'(0) = fix(0) > 

which implies that y(t) > for all t € (0, e) for some small e > 0. The fact that this holds at all times follows 
from the same reasoning as before based on the fact that both functions are bounded. H 
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The next lemma, which also follows from Lemma [21 is our main tool to prove Theorems [T]|H It lists some of 
the invariant sets of the system. 

Lemma 4 The following sets are positive invariant for the system JU-JH]). 



fie 


■= {{x,y) e Qq 


: x > 9 and 


2/>0} 




■= {(z,y) G Q 


: x > 1 and 


y> 1} 


fio,e 


■= {(x,y) g fi 


: < .t < 6» 


and <y <9} 


fio,i 


■= {(x,y) g I2 


: < x < 1 


and 9 < y < 1} 


^x<y 


:= {(a;,?/) G fi Q 


: x < y} 




^x>y 


:= {(x,y) G 


: x > y} 




^x—y 


:= {(a;,y) £ I2 


: x = y}. 





Moreover, the dynamics along the invariant set fi x=y are described by 

1. If xq = i/o€ (0,0) then x(t) = y(t) — > as t — > oo. 

2. If xq = yo G (0, oo) f/ien x(£) = y(t) —> 1 as i — > oo. 

Proof First, we assume that the initial condition [xo,yo) G fig and introduce 

u(t) = x(t) - 9 and v(t) = y(t) - 9. 

Then, the system ((I])-© can be rewritten as 

u = x = r u (u + 9)(1 — 9 — u) — fi (v — u) (8) 
v = y = r v (v + 9)(1 — 9 — v) — fj,(u — v) (9) 

with initial condition (uq, vq) G J?o- Now, the arguments of the proof of Lemma [2] imply that fig is positive 
invariant for (JS])-©- Moreover, 

(it(i), v(t)) e fi Q if and only if (x(t), y{t)) G fig 

so the set fig is positive invariant for the system (@|-([5]). The fact that fl\ is positive invariant follows from 
the same argument but applied to 

u(t) = x(t) - 1 and v(t) = y(t) - 1. 

To prove the positive invariance of fio t g we first observe that Lemma [2] implies that any trajectory starting 
from a point in the square fig,g cannot exit the square crossing its left of bottom side. Moreover, the same 
arguments as in the proof of Lemma [2] imply that it cannot exit the square crossing its right or top sides 
either because of the following three properties. 

1. The upper right corner (9,9) is a fixed point of the system (JU)-©. 

2. If (x(t),y(t)) G fi Q .g with x{t) = 9 then 

x'(t) = »(y(t)-x(t)) = n{y{t)-6) < 0. 

3. If (x(t),y(t)) G fi M with y(t) = 9 then 

y'(t) = (x(x(t)-y(t)) = fi(x(t)-9) < 0. 
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This proves that J?o,e is positive invariant. The fact that the square is also positive invariant follows 
from the same argument, looking at the derivatives along each side and using that the four corners are 
equilibriums. To prove the positive invariance of the last three sets, we introduce the new functions 

x(t) + y{t) x(t) - y(t) 
u{t) = and v(t) = - . 

A straightforward calculation shows that 

ii = r u(u — 0)(1 -u) + rv 2 {l + 6-3u) (10) 
v = r v (-3u 2 + 2(1 + 8)u~v 2 - 9 - 2/i/r). (11) 

From ()ll[l. we see that v = is an invariant manifold of v, i.e, 

v(0)=0 implies v(t)=0 for all t > 0, 

from which it follows that the set fl x=y is positive invariant for the original system (H])-©. In particular, if 
Vq = then (fT0|) reduces to 

ii = ru (u — — u). 
Therefore, by applying Lemma [TJ we can conclude that 

1. If Mo £ (0, 6) then u(t) — > as t — >• oo, 

2. If ito G (0, oo) then u(t) — > 1 as t — >• oo, 

which, in view of the definition of it and w, and the fact that fl x=v is positive invariant, is equivalent to the 
last two statements of Lemma [4] Finally, for any initial condition xq > ya, Lemma [2] implies that u(t) and 
v(t) are both bounded uniformly in time so, using equation (|11[) and the same argument as in the proof of 
Lemma [31 we can deduce that 

x(t)-y(t) = 2v(t) > 2v exp(-Kt) = (x - y ) exp(-Kt) > 

for alH > and some constant K < oo. This proves that fi x >y is positive invariant. By symmetry, the same 
holds for the set f2 x <y M 

With Lemmas [2HH in hands, we are now ready to prove the main results for the deterministic two-patch 
model described by the system (H])-©. 



Proof of Theorem [T] 



By Poincare-Bcndixson Theorem, the omega limit set of the system ((¥])- {5) is either a fixed point or a limit 
cycle. If the inequality © holds, we can use Dulac's criterion to exclude the existence of a limit cycle. Let 
c G [0, 3) and define the scalar function p c (x, y) = (xy)~ c on R^_. Then, 

d 

— \{rx (x - 0)(1 -x)+n(y- x))p(x, y)\ 

d 

+ l(ry (y - - y) +n(x- y))p(x, y)\ 

= (xy)- c [r(c - 3)(.t 2 + y 2 ) + r(2 - c)(l + 6)(x + y) 



+ 2r6(c -l)-2/j, + cfi(2 - xy- 1 - yx' 1 )] 



< (xy)~ 



r(c — 3) ( x + 



(2-c)(l + , 



2(c-3) 



r(2-c) 2 (l + 6) 2 



4(c-3) 



+ r(c-3) y + 



(2-c)(l + . 
2(c-3) 



r(2 - c) 2 (l + i 
4(c-3) 



+ 2rd{c- 1) - 2^ 
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In particular, if © holds then the equation above is strictly negative for any [x, y) G Q \ E . Therefore, by 
Dulac's criterion, the system has no limit cycle, i.e., any trajectory of (HJ)-© starting with a nonncgative 
initial condition converges to a fixed point. ■ 



Proof of Theorem [3] 

Parts 1 and 2 of Theorem [2] about the local stability of the three symmetric equilibriums follow from the 
analysis of the Jacobian matrices. For each of the three equilibriums, we have 

Eq - The Jacobian matrix associated with this equilibrium is 

f-rO-fi n \ 
Jo = (12) 
\ M -rO -fx J 

with eigenvalues Ai = —rO and A2 = —rO — 2/i associated with (1, 1) and ( — 1, 1) as their eigenvectors, 
respectively. We can easily conclude that the trivial boundary equilibrium Eo is locally stable since both 
eigenvalues of (fT2"]) are negative. 

Eg - The Jacobian matrix associated with this equilibrium is 

fr9(l-6)-fx fx \ 

Je = (13) 
V H r6(l-6)J 

with eigenvalues Ai = r6(l — 0) and A2 = r9(l — 0) — 2/i associated with (1,1) and (—1,1) as their 
eigenvectors, respectively. We can easily conclude that the equilibrium Eg is always unstable on the 
invariant set f2 x=y . Moreover, if 2/i > r9(l — 9) then Eg is a saddle, while if 2/i < r9(l — 9) then Eg is a 
source. 

Ei - The Jacobian matrix associated with this equilibrium is 

Ji = [ (14) 

with two negative eigenvalues Ai = — r(l — 9) < and A2 = — r(l — 6) — 2/i < since 6 < 1. Therefore, 
the equilibrium E\ is also locally stable. 

To prove the third part of the theorem, we first define the function u(t) = x(t) + y(t). Then 

u = x + y = rx(x — 9)(l — x) + ry (y — 9)(l — y). 

To prove that 

f2 Oi a\{(0,0)} C B Q (15) 

we first assume that 

(xo,yo) € flo,6 \ {E ,E e }. 
Since the set f2o,e is positive invariant, we have u'(t) < for all t > 0. Using in addition that 

u{t) = if and only if u(t) = 0, 

we can conclude that u(t) converges to zero. Recalling the definition of u and invoking again the positive 
invariance of i^o,e, we can deduce that x(t) and y(t) converge to zero so (|15[) holds. To prove that 

n e \{{0,0)} c b 1 (16) 

we now assume that 

(x ,y ) e i7 e \ {Eg, E x }. 

Then, we have the following alternative. 
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1. (xo,yo) £ ^0,1 \ {Eg, Ei}. Since fie,i is positive invariant, we may use the same argument as before to 
see that the derivative of u is nonncgative and the system converges to the equilibrium point E\ . 

2. (xo, uq) £ Q\ \ {E±}. Repeating again the same argument but with the positive invariant set fl\ implies 
that the system converges to E\ . 

3. (xo, yo) £ fig \ {fte,i U We may assume that xq < yo without loss of generality since the system is 
symmetric. Then, using the positive invariance of the set fl x<y we have x(t) < y(t) for all t > so 

x = rx (x — 9)(1 — x) + fi (y — x) > if x < 1 

V = ry (y - - y) + /i (x - y) < ify>l. 

This indicates that the trajectory starting at (xq, yo) can only exit the infinite rectangle [9, 1] x [1, oo) by 
crossing its bottom or right side. Therefore, we have the following three possibilities. 

a. No exit: (x(t), y(t)) £ Qg t i\Jfli for all t > 0. In this case, the sign of the derivatives implies convergence 
to Ei. 

b. Bottom side: (x(t),y(t)) £ for some time t > 0. In this case, point 1 above implies convergence 
to the equilibrium point E\ . 

b. Right side: (x(t),y(t)) £ i?i for some time t > 0. In this case, point 2 above implies convergence to 
the equilibrium point E\ . 

Combining 1-3 above implies (|16|) . Now assume that < 1/2. Defining 

u(t) = X{t)+ 2 y(t) and „(*) = X{t) ~ yit) 

recall that the system ((5J can be rewritten as 

ii = ru(u-6)(l-u) + rv 2 (l + 9-3u) 
v = r v {-'iu 2 + 2(1 + 6)u-v 2 - 9 - 2/x/r). 

To prove that 

B C {(x,y) £ n :x + y < 29} (17) 

it suffices to prove that 

xo + Vo > 29 implies x(t) + y(t) > 29 for all t > 0. (18) 

Assume by contradiction that (|18|) is not satisfied. Then, there exists an initial condition with xo + yo > 20 
and a time T > such that 

x(0) + y(0) = x + yo > 29 and x{T) + y(T) < 29. 

By continuity of the trajectories, the intermediate value theorem implies the existence of a time t < T such 
that u(t) = 29 therefore 

S := sup{t < T : u{t) = 29} 

is well defined and u(t) < 29 for all t £ (S,T]. To prove that this leads to a contradiction, we consider the 
following two cases. 

1. If xo 7^ yo, the invariance of fl x <y and f2 x>y implies that x(t) ^ y(t) at any time, from which it follows 
that 

u'(S) = rv 2 (S){l + 9 -39) = (r/4) (x(S) - y(S)) 2 (l - 29) > 0. 
In particular, there exists e > such that u(t) > 29 for all t £ (S,S + e), which contradicts the existence 
of time S. 

2. If xo = yo the result directly follows from the fact that fig n fi x =y is positive invariant, as it is the 
intersection of two invariant sets. 

Combining 1 and 2 above yields (fP7|) . The proof of the last inclusion in Theorem [2] follows from similar 
arguments. ■ 
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Proof of Theorem [3] 

Define as previously the functions 



Using (TTTj) above, we obtain 



u(t) = x ®+y® and V (t) = AzM. 



1 + ey 2 / r(0 2 - 6>+ i) 



r v —3 u — /i — v 



r 



,2 



Assume first that xq > yo- Using the fact that the set fl x > v is positive invariant by Lemma 21 we deduce 
that x(t) > y(t) at all positive times t. In particular, recalling the definition of v, and using the expression 
of the derivative v above and the fact that ([7]) holds, we obtain that 

v(0) > implies v(t) > and v'(t) < for all t > 0. 

Since v = if and only if v = 0, we deduce that v(t) converges to 0. By symmetry, the same can be proved of 
the system starting with any initial conditions such that xq < yo. Since £2 x=y is positive invariant, we have 
the same conclusion when the initial condition satisfies xq — yo, which can also be seen from the expression 
of the derivative v. Therefore, if ([7]) holds then 

for all (x ,yo) G Ki, lim v(t) = 0, 

t— > oo 

so for any (xo, yo) G and any e > 0, there exists k > such that 

\x{t) - y(t)\ < e for all t > k. 

It follows that any trajectory of the system converges to one of the symmetric equilibriums Eo, Eg or E\. 
Now, observing that 

r{6 2 -e + l) r0(l-0) r(20-l) 2 



6 2 6 

we obtain that 

r(0 2 -0 + l) r0(l- 



> 0, 



A 4 > a - 



6 " 2 

In view of the expression of the Jacobian matrix (|13[) . this implies that the equilibrium (0, 0) is a saddle with 
unstable manifold 

Qx=y \ {Eo, Eg, Ei}. 

In particular, it follows from Hartman-Grobman Theorem and the second part of Lemma 2] that there are 
only two attractors: Eq and E\. Hence, the system starting from any initial condition not belonging to the 
manifold Sg converges to either E or E\. To conclude the proof, it suffices to observe that, by (TTS)) in the 
proof of Theorem [H if < 1/2 then the set 

{(x,y) eQ \Se :x + y>26} 

is positive invariant so the system starting from any initial condition in this set converges to E\, the only 
attractor in this invariant set. Similarly, the last statement follows from the fact that, if > 1/2 then the set 

{(x,y) e Q\S :x + y<29} 

is positive invariant and contains only one attractor: Eq. M 
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Fig. 5 Schematic presentations of nullclines of the system with black dots representing locally stable equilibriums. 

Proof of Theorem [4] 

We first prove that all the equilibriums of the system (jl])-© belong to the unit square [0, 1] x [0,1]. Since 
the system is symmetric and, by Lemma01 the omega limit set of any initial condition (xo,yo) £ Qo,6 U fig 
belongs to the unit square (either Eq, Eg or E%), it suffices to focus on the case 



Since xq = ya only happens when starting from (0,0), to avoid trivialities, we shall assume in addition that 
xq < 2/o- Then, applying Lemma [4[ we obtain that x(t) < y(t) for all i > which, together with ([5]), implies 
that 



Excluding the trivial case when the initial condition belongs to the stable manifold of (9,9), in which case 
its omega limit set reduces to (9,9), we have the following alternative. 

1. If (x(t),y(t)) e fig for some t > then, by the second part of Theorem [21 the omega limit set of the 
initial condition is E\. 

2. If (x(t),y(t)) ^ fig for all t > then ([TO)) and the fact that f2 x<y is positive invariant imply that 
x(T) < y(T) < 1 for some T > 0. Using as previously the continuity of the trajectories and the fact that 



allows to invoke the intermediate value theorem and prove by contradiction that x(t) < y(t) < 1 at any 
time t > T. 

This establishes the first part of Theorem [U The second part follows directly from the proof of Theorem [3] 
To prove the third part, we first observe that the equation of the nullclinc y = f(x) can be rewritten as 



In particular, if (1 — 9) 2 > A^/r then the nullcline intersects the x-axis at the three points with coordinates 
(0,0), (xi,0) and (x2,0) where 



(x , Vo) 6 {(x, y)en :0<x<9<y}. 



V = ry(y -6)(l-y) + n(x-y) < ify>l. 



(19) 



y = V (x - y) < if y = 1 






Finally, a phase-plane analysis based on Figure [5] shows that 



26 



1. If x\ < M < X2 then the system has five fixed points with only two locally stable equilibriums: Eq and 
El 

2. If M > 1 then the system achieves maximum number of equilibriums which is nine, with only four locally 
stable equilibriums. 

3. If the system has less than nine equilibriums, then it has only two local stable equilibriums: Eq and E\. 
Note also that M can be computed explicitly: M = f(x*) where 

x* = i ((1 + 9) - - 9 + 1 - 3///r) 
is the smallest root of the polynomial h'(x). ■ 

Proof of Theorem [5] 

The first step is to prove that the set of upper configurations is closed under the dynamics. We observe that, 
condition on the event that the configuration is an upper configuration, only expansions and migrations can 
occur. Furthermore, migration events can only result in an increase of the lowest density and a decrease of 
the highest density, i.e., if a migration event occurs at time t and the configuration at time t — At is an upper 
configuration then 

mm(X t - At, Y t _ At ) < mm(X t ,Y t ) 

< max(X u Y t ) < max(X i _ 44 ,F t _ /it ). 

It follows that the set of upper configurations (and similarly the set of lower configurations) is closed under 
the dynamics, i.e., 

p((x t ,Y t ) e n+ | {x ,Y Q ) g n+) 

= P((x t ,Y t ) g n- | (x ,y ) en-) = l 

for all times t. Since, starting from an upper configuration, the system jumps to (1, 1) whenever two expansion 
events at X and Y occur consecutively (they are not separated by a migration event), we deduce that the 
stopping time r + is almost surely finite. The same holds for the stopping time t~ when starting from a lower 
configuration. Hence, 

P(r = r+ < oo | (X ,Y ) G Q+) 

= P{t =t~ < oo | (X ,Y ) G Q-) = 1. 

To compute the expected value of the time to fixation, we now construct the stochastic process graphically 
from a collection of Poisson processes, relying on an idea of Harris (1972). Two Poisson processes, each with 
parameter r, are attached to each of the patches X and Y, and an additional Poisson process with parameter 
one is attached to the edge connecting the patches. All three processes are independent. Let 

r x = {T* : n > 1}, TV = {T„ Y : n > 1}, r e = {T n e : n > 1} 

denote these Poisson processes. At any time of the process rx the population size at patch X jumps to 
cither or 1 depending on whether it is smaller or larger than 9 by this time, respectively. The evolution 
at patch Y is defined similarly but using the Poisson process TV- At each time in r e , a fraction fx of the 
population at each patch is displaced to the other patch. To compute the expected value, we let t > and 
introduce the stopping times 

T z = min {r z n (t, oo)} for Z = X, Y, e. 

Then, P (max(Tx, Ty) < T e ) is the probability that two consecutive migration events are separated by at 
least one extinction-expansion event at patch X and one extinction-expansion event at patch Y. To compute 
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this probability, we first observe that Tx and T Y are independent exponentially distributed random variables 
with parameter r, from which it follows that 

P (max(Tx, T Y ) < u) = P (T x < u, T Y < u) = (f - cxp(-ru)) 2 . 

Since T e is exponentially distributed with parameter 1, 

P(max(T x ,T Y ) <T e ) = J J e~ v — [{1 - exp(-ru)) 2 J dv du 

2r 2r 2r 2 

■= Ps- 



r+l 2r+l (r+l)(2r+l) 

Hence, the last time a migration event occurs before fixation is equal in distribution to T|_ 1 where the 
random variable J is a geometrically distributed with parameter p s from which we deduce that 

E [t+ I 9 < X , Y < 1] = E [T e ] x E [J - 1] + E [max(T X) T Y )] 
(r + l)(2r + l) 



2r 



2 



/■CX) 

1+ / P(max(T x ,T Y ) > w) 



3r + ! r , , nn 2j 6r + l 

= + X 1 - (! - CX P(~™)) du = -^S-- 

The same holds for the stopping time r~ when starting the process from a lower configuration. This completes 
the proof of Theorem [5] ■ 



Proof of Theorem [S] 

We first prove that P (T < oo) = P (r < oo) = 1. Let e > small. Then, for almost all realizations of the 
process, there exists an increasing sequence of random times T\ < • • • < Tj < • • • such that 

lim T t = oo and \X Ti + Y Tl - 29\ > e for all i > 1. 

Moreover, there exists K < oo that does not depend on i such that, if after T,; a sequence of ii" migration 
events occur before any expansion or extinction events then the system hits either an upper configuration or 
a lower configuration. Since K is finite, such an event has a strictly positive probability, so the Borel-Cantelli 
Lemma implies that the process hits either an upper configuration or a lower configuration after a random 
time which is almost surely finite: P (T < oo) = 1. Theorem [5] then implies that 

P (t < oo) = P (t+ < oo) + P(t~ < oo) 

> P (t+ < oo I (X , Y ) € 0+) P (T+ < oo) 

+ P (t- < oo | (X , Y ) € Q-) P (T~ < oo) 
= P (T+ < oo) + P (T~ < oo) = P (T < oo) = 1. 

To estimate the expected value of T, we observe that the transition rates of the process indicate that if at 
time t exactly n migration events but neither expansion nor extinction events have occurred then 



(X t ,Y t ) = f n (0 : l) where f(a, b) = (1 - /i) (a, b) + /i (6, a), 
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(0, 1) Ui u i+ \ (0, 1) 




Vn 2 

Fig. 6 Schematic representation of the stochastic process (Ut, Vt). 



so that X t < u n and Y t > v n where u n and v„ are defined recursively by 

u n+ i = (1 — fi) u n + /i with u a = 0, 
v n+ i = (1 - n) v n with v a = 1. 

A straightforward calculation shows that 

n— 1 n— 1 

u„ = (u fc+ i - u fe ) = (1 - ^) fc (in - u ) = 1 - (1 - /i)" 

fc=0 fc=0 

and u n = 1 — u n = (1 — fj) n , therefore 

u n > 8 if and only if n > n\ := [ln(l — 0)/ln(l — /i)J 
u„ < 6 if and only if n > n-i := [ln(^)/ ln(l — /i)J 

where [-J is for the integer part. Now, let {(Ut,Vt)}t be the Markov process with state space 

E = {(ui,Vj) >0} 

and transition rates 

P {(Ut + AuV t+At ) - (0, Vj ) | (U t ,V t ) = (ui,Vj)) - rAt + o(At) 
P {{Ut+AuVt+At) = (Ui, 1) | {U t , V t ) = (ui, Vj )) = rh + o(h) 
P ' ((Ut+AuVt+At) = | (U t ,V t ) = (iii,Vj)) = At + o{At) 

and starting at (Uq, Vq) = (0, 1). We call W-, N-, and S'iJ-jumps, the jumps described by the three transition 
rates above, respectively, and refer the reader to the left-hand side of Figure[6]for an illustration of the process. 
By construction of the sequences {u n ) n and (v n ) n , we have 
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P (X t > a | T > t) < P(U t >a) 
P(Y t >a\T>t) > P(V t >a) 

for all a G [0, 1], i.e., before the process hits an upper or a lower configuration, X t is stochastically smaller 
than Ut while Y t is stochastically larger than Vt. This implies that E [T] > E [T*] where 

T* = inf {t > : U t > u ni or V t < v n2 }. 

Ei = {{ui, Vj) : i,j < uq} and E2 = {(ui, Vj) : i < m and j < ni\. 
Then, T* is the first time (Ut, Vt) exits the set E2, i.e., 

T* = mi{t>0:(U t ,Vt)iE 2 }. 

So, to bound E [T*] from below, it suffices to prove that (Ut, Vt) £ E2 for an arbitrarily long time. The idea 
is to prove that, when starting from the smaller rectangle E\, the process stays in E2 and comes back to 
Ei after no jumps with probability close to 1. Using in addition the Markov property, we obtain that the 
number of jumps required to exit E2 is stochastically larger than no times a geometric random variable with 
small success probability. To make this argument precise, we let (U n , V n ) denote the embedded discrete-time 
Markov chain associated with the process (Ut, Vt). To count the number of steps needed to exit the rectangle 
E2, we define a sequence of Bernoulli random variables {Z^ : k > 1} associated to (U n , V n ) by setting 

Zk = if there is at least one iV-jump and one VF-jump 

between time (k — 1)uq + 1 and time kno 

= 1 if there is no TV-jump or no W^-jump 

between time (k — 1)uq + 1 and time kno 

Since (U n ,V n ) is a discrete-time Markov chain, the random variables Zk are independent Bernoulli random 
variables, and a straightforward calculation shows that the success probability is given by 

( r \ no fl + r\ n " 

P(Z k = l)<2[l =2^^— . 

Moreover, since no = (1/2) min(ni,n2), we have that 

(Ukn , Vkn ) S Ei and Z k+ i = implies that 

(Un, V n ) 6 E 2 for all kn < n < (k + l)n and (U (k+ i )na , V( fc+ i)„ ) G E x . 

See the right-hand side of Figure |6l This indicates that 

Zi = Z 2 = ■ ■ ■ = Z k = => (U n , V n ) G E 2 for all n < kn . 

Finally, using that (Ut, Vt) jumps at rate 1 + 2r and that inf {k : Z k = 1} is stochastically larger than a 
geometric random variable Z with success probability P (Z k = 1) we can conclude that 

E [ T] > E[r] > _J!2_e[Z] = -52- fl±^V°. 
1 i ~ 1 1 ~ 1 + 2r 11 2 + 4r \l + r J 

This completes the proof. ■ 
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Proof of Theorem [7] 



First, we observe that the process Ut introduced in the proof of Theorem [5] is stochastically larger than X t 
so to prove the first inequality it suffices to establish its analog for the expected value ¥. v (Ut) where n is the 
stationary distribution of the stochastic process Ut ■ Note that the infinitesimal matrix of the Markov process 
Ut expressed in the basis (mq, u\, v,2, ■ • •) is given by 



Q 



f-l 1 o o •■•^ 

r -(r + 1) 1 ■•■ 

r -(r+1) 1 

r '■• 



By solving tt ■ Q = 0, we find that 



r + 1 V r + 1 



r + 1 



r + 1 



This implies that 



n=0 



r + 1 



71+1 



I 1 

E„ (X t ) < E x (U t ) = r V u n [ — 

Z — i \r + 1 

oo 

= r ^ (l-(l-M) 

71=0 

00 / 1 

= — y f— 1 - — y 

r + l^Vr + l/ r+l^ 



71+1 



1-M 
r + 1 



r + /i 



The proof of the second inequality is similar. 



Proof of Theorem [8] 

We first observe that the processes (X t ,Y t ) and (X t ,Yt) can be constructed on the same probability space 
starting from the same initial configuration in such a way that X t = X t and Y t = Y t until the hitting time 
T, which we assume from now on. Let Tq = and, for all % > 1, let Tj denote the time of the ith jump of 
the process £ t := X t + Y t . Since migration events do not change the value of £t, time T, corresponds to the 
time of an extinction event at X or an expansion event at Y, therefore we have 

T = T + if and only if there exists i > 

such that T 6 (Ti,T i+ i) and > 26». 

T = T~ if and only if there exists i > 

such that T 6 (Ti,T i+1 ) and < 20. 

Let e > small such that 1 — e > 29, and consider the events 



D i,n = i X Ti = and Y Ti £29- [ne, (n + l)e)} 
D tn = {^i - 1 and X Ti e [ne, (n + l)e)}. 
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First, since 1 — (n + l)e > 29 — ne, migration events between Tj and T, + i displace less individuals on the 
event D~ n than on Df n so 

P(T G (r <} r 4+1 ) | Dr n ) < P (T G (T,,T l+1 ) I £>+J. 

Second, note that ti n = (1 — /i)" < 20 if and only if we have 

n > m := [hi(29)/ ln(l - #)J . 

In particular, if is the time of an extinction event at X then ij^ < 20 only if at least mo migration events 
have occurred since the last expansion event at patch Y . This implies that 

P (Y Tz < 20) < 



1 + r 



Since by symmetry the random variables X t and 1 — Y t are identically distributed, and 29 — X t is stochastically 
smaller than Y t , we deduce that 



P(DrJ < P(Y Ti <26) P(Df n ) < I — j P(D+ n ). 

Finally, observing that 

{T = T-} = (J |J {T G (T h T i+1 )} n D~. 

i=0 n=0 

oo le- 1 } 

and {T = T+} D |J |J {Te (T„T, +1 )} n 

i=0 n=0 



we can conclude that 



oo L^J 

P(T = T-) = J2 E F ( T£ ( T i' T i+l) I D n,i) P ( D n,i) 
i=0 n=0 

oo [e- 1 ] 

<EE P{TG(T h T i+1 ) | £>+,) P (£»+,) (P{D- ti )/P (D+.)) 

i=0 n=0 

/ i \ mo oo L £ J 

* (lT7 E E P{T£{T h T i+1 );D+J 

i=0 n=0 

m 



1 + r 



< (t— ) P{T = T+). 



This completes the proof. 
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